import pandas as pd
import numpy as np
import os
import shap
import xgboost as xgb
from sklearn.linear_model import Ridge
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import precision_score, recall_score, accuracy_score, roc_curve, \
roc_auc_score, precision_recall_curve, confusion_matrix, r2_score, mean_absolute_error, \
mean_squared_error
import pandas_profiling as pp
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:60% !important; }</style>"))
Read the data
!dir data
data_dir = 'data'
df = pd.read_csv(os.path.join(data_dir, 'MEPS_data_preprocessed.csv'))
df.reset_index(drop=True, inplace=True)
df_raw = df.copy()
df.head()
Change categorical variables to strings (even though they're represented as numbers)
cats = ['REGION','MARRY31X','EDRECODE','FTSTU31X','ACTDTY31','HONRDC31',
'RTHLTH31','MNHLTH31','HIBPDX','CHDDX','ANGIDX','MIDX','OHRTDX','STRKDX',
'EMPHDX','CHBRON31','CHOLDX','CANCERDX','DIABDX','JTPAIN31','ARTHDX',
'ARTHTYPE','ASTHDX','ADHDADDX','PREGNT31','WLKLIM31','ACTLIM31','SOCLIM31',
'COGLIM31','DFHEAR42','DFSEE42','ADSMOK42','PHQ242','EMPST31','POVCAT15','INSCOV15']
df = df.copy()
for col in cats:
df[col] = df[col].astype(str)
# # Show summary report
# df.profile_report(title='Summary Report')
Let's have a closer look at the predicted variable (HEALTHEXP)
df['HEALTHEXP'].hist(bins=100)
df['HEALTHEXP'].describe(percentiles=np.linspace(0,1,21))
health_exp = df['HEALTHEXP'].values
plt.plot(np.linspace(0,100,21), np.percentile(health_exp, np.linspace(0,100,21)),'o');
plt.axhline(y=np.mean(health_exp), color='red', label = 'mean')
plt.xlabel('Percentile');
plt.ylabel('Healthexp');
plt.legend();
plt.plot(np.linspace(0,100,21), np.percentile(health_exp, np.linspace(0,100,21)),'o');
plt.axhline(y=np.mean(health_exp), color='red', label = 'mean')
plt.xlabel('Percentile');
plt.ylabel('Healthexp');
plt.yscale('log')
plt.legend();
plt.figure(figsize=(6,6))
df[['HEALTHEXP']].boxplot();
plt.figure(figsize=(6,6))
plt.boxplot(df['HEALTHEXP']);
plt.title('Health Exp')
plt.yscale('log')
Notes on predicted variable:
cut_off = 100_000
df = df[df['HEALTHEXP']<cut_off]
Prepare features
# Drop panel number (not meant to be predictive) and sample weights
df.drop(columns = ['PANEL', 'PERSONWT'], inplace=True)
df.head()
y = df.pop('HEALTHEXP')
One-hot encoding for all the categorical features
cat_cols_one_hot = []
for col in cats:
var_one_hot = pd.get_dummies(df[col],prefix=col, drop_first = True)
df = df.drop(columns=[col])
df = pd.concat([df, var_one_hot], axis=1)
for cl in var_one_hot.columns:
cat_cols_one_hot.append(cl)
df.head()
More feature engineering in next HW ...
def train_print_save_models_reg(mods, dfTrain, dfTest, yTrain, yTest):
trained_models = {}
for model_name, mod_object in mods.items():
arTrain = np.asarray(dfTrain)
arTest = np.asarray(dfTest)
md = mod_object.fit(arTrain,yTrain)
y_pred_test = md.predict(arTest)
y_pred_test[y_pred_test<0]=0
y_pred_train = md.predict(arTrain)
mae_train = mean_absolute_error(yTrain, y_pred_train)
mae_test = mean_absolute_error(yTest, y_pred_test)
rmse_train = np.sqrt(mean_squared_error(yTrain, y_pred_train))
rmse_test = np.sqrt(mean_squared_error(yTest, y_pred_test))
r2_train = r2_score(yTrain, y_pred_train)
r2_test = r2_score(yTest, y_pred_test)
print(model_name)
print('Training R^2: {}, MAE:{}, RMSE:{} '.format(r2_train, mae_train, rmse_train))
print('Test R^2:{}, MAE:{}, RMSE:{}'.format(r2_test, mae_test, rmse_test) )
trained_models[model_name] = md
return trained_models
mods = {
'XGB':xgb.XGBRegressor(max_depth=3),
'LinReg':LinearRegression(),
'Ridge_Reg':Ridge(alpha=.5)
}
dfTrain, dfTest, yTrain, yTest = train_test_split(df, y, random_state=42)
len(dfTrain), len(dfTest)
models = train_print_save_models_reg(mods, dfTrain, dfTest, yTrain, yTest)
import tensorflow as tf
from keras.models import Sequential
from keras.optimizers import RMSprop
from keras.layers import Dense, Dropout, Embedding, Reshape
from keras.backend import clear_session
from keras.optimizers import Adam
from keras import backend as K
from keras.models import Model
from keras.layers import Dense
clear_session()
def r2(y_true, y_pred):
SS_res = K.sum(K.square( y_true-y_pred ))
SS_tot = K.sum(K.square( y_true - K.mean(y_true) ) )
return ( 1 - SS_res/(SS_tot + K.epsilon()) )
model = Sequential([
Dense(64, activation=tf.nn.relu, input_shape=[dfTrain.shape[1]]),
Dropout(0.3),
Dense(128, activation='relu'),
Dense(64, activation='relu'),
Dense(1)
])
optimizer = Adam(1e-4)
model.compile(loss='mae',
optimizer = optimizer,
metrics=[r2])
model.fit(x=dfTrain,
y=yTrain, epochs=1000,
batch_size =1024,
validation_data=(dfTest,yTest),
verbose=2,
validation_freq=1, use_multiprocessing=True)
Simple ANN
y_pred = model.predict(dfTest)
r2_score(yTest, y_pred)
mae = mean_absolute_error(yTest, y_pred)
rmse = np.sqrt(mean_squared_error(yTest, y_pred))
print('MAE: {}, RMSE: {}'.format(mae, rmse))
models['ANN'] = model
arTrain = np.asarray(dfTrain)
arTest = np.asarray(dfTest)
model1 = models['XGB']
model1
obs_num_1 = 10
arobs = np.reshape(arTest[obs_num_1], (1,-1))
dfTest.iloc[obs_num_1,:]
#wartość rzeczywista
yTest.iloc[obs_num_1]
#predykcja
model1.predict(arobs)[0]
import lime
from lime.lime_tabular import LimeTabularExplainer as LTE
cat_features = np.asarray([i for (i,col) in enumerate(dfTrain.columns) if col in cat_cols_one_hot])
feat_names = np.asarray(list(dfTrain.columns))
dfTrain.columns[cat_features]
explainer = LTE(arTrain, mode='regression', feature_names=feat_names,
class_names=['asd'], categorical_features=cat_features,
verbose=True)
exp = explainer.explain_instance(arTest[10], model1.predict, num_features=6)
exp.show_in_notebook()
Opis zmiennych wybranych przez explainer i możliwe wartości:
Uwagi:
Model jako najistotniejszą cechę wybrał fakt bycia przez pacjenta w ciąży (w okresie, którego dotyczy zebrany wywiad). Wydaje się to mieć sens, jeśli kobieta była w ciąży to na pewno korzystała z usług medycznych - cały okres okołociążowy obiftuje w wizyty i badania lekarskie. Inna sprawa, że ta wybrana obserwacja dotyczy mężczyzny.
Druga wskazana cecha to fakt, że wskaźnik określający ogólny stan zdrowia fizycznego nie mógł zostać wyliczony (lub jest bardzo niski?). Być może zmienna ta wymaga modyfikacji.
Następnie model wskazuje na fakt niezdiagnozowania u pacjenta cukrzycy => wydatki na OM mniejsze.
Wyjaśnienie mówi też, że brak ograniczeń w wykonywaniu prac zarobkowych/domowych ujemnie wpływa na wydatki na opiekę medyczną.
Fakt niezdiagnozowania nowotworu ujemnie wpływa na przewidywane wydatki na OM.
Trudno jest zinterpretować wskazania modelu w przypadku samooceny stanu zdrowia (RTHLTH31). Wyjaśnienie mówi, że jeśli wartość ta jest różna od 4, to wydatki są generalnie mniejsze.
df_raw['RTHLTH31'].value_counts()
dfTest.iloc[obs_num_1,:]
df_raw['PHQ242'].hist(bins=2*7);
plt.plot()
several_obs = [23,32,42]
dfTest.iloc[several_obs,:]
for obs_id in several_obs:
explainer.explain_instance(arTest[obs_id], model1.predict, num_features=6).show_in_notebook()
Zmienne, które nie pojawiły się wcześniej:
Widzimy, że wśród tych kilku wybranych obserwacji pojawia się zestaw powtarzających się cech i ich wpływ na model (jakościowo) jest taki sam. Najistotniejsza pozostaje zmienna dotycząca ciąży. Wszystkie zmienne, które odpowiadają za fakt zdiagnozowania/niezdiagnozowania konkretnych chorób (kończące się na DX) wydają się sensownie wpływać na wyjścia modelu (jeśli ktoś ma/miał zdiagnozowaną daną chorobę to wydatki rosną). Dodatkowo pojawiają się zmienne związane z samopoczuciem i oceną stanu zdrowia, których interpretacja znów jest utrudniona.
Jeden z modeli wskazał, że fakt iż osoba na pewno nie jest doktorem/magistrem (EDRECODE_16=0) powoduje, że wydatki na OM są mniejsze.
Uwaga ogólna - kodowanie jedynkowe utrudnia odczyt takich wykresów.
model2 = models['ANN']
model2
ar = np.asarray(dfTest)
model2.predict(arobs)
exp1 = explainer.explain_instance(arTest[10], model1.predict, num_features=6)
exp2 = explainer.explain_instance(arTest[10], model2.predict, num_features=6)
exp1.show_in_notebook()
exp2.show_in_notebook()
Zmienne, które nie pojawiły się wcześniej:
df_raw['RACE3'].value_counts()
for obs_id in several_obs:
explainer.explain_instance(arTest[obs_id], model2.predict, num_features=6).show_in_notebook()
Porównanie modeli dla obserwacji nr 10 i kilku losowych: